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Abstract 

We present a family of p-enrichment schemes. These schemes may be separated into two basic 
classes: the first, called fixed tolerance schemes, rely on setting global scalar tolerances on the local 
regularity of the solution, and the second, called dioristic schemes, rely on time-evolving bounds on 
the local variation in the solution. Each class of p-enrichment scheme is further divided into two basic 
types. The first type (the Type I schemes) enrich along lines of maximal variation, striving to enhance 
stable solutions in "areas of highest interest." The second type (the Type II schemes) enrich along 
lines of maximal regularity in order to maximize the stability of the enrichment process. Each of these 
schemes are tested over a pair of model problems arising in coastal hydrology. The first is a contam- 
inant transport model, which addresses a declinature problem for a contaminant plume with respect 
to a bay inlet setting. The second is a multicomponent chemically reactive flow model of estuary 
eutrophication arising in the Gulf of Mexico. 
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namic p-enrichment, flow reactors, advective transport, shallow water equations, contaminant trans- 
port, contaminant declinature, estuary eutrophication. 
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§1 Introduction 

The topic of p-enrichment holds substantial allure in computational engineering and numerics due to 
the promise of being able to locally resolve solution structure without paying the full cost of raising 
the global order of the solution a complete integral step. Though a good number of important results 
exist regarding p-enrichment, the majority of these results deal with p-adaptation as the fundamentally 
coupled component of an /ip-adaptive scheme, wherein the numerical behavior of the p-enrichment is not 
isolated from the behavior of the /i-adaptivity. Moreover, of the work that does address p-enrichment 
as its own topic, there is often the strong presence of assumptions restricting the resulting behaviors 
to highly idealized linear systems of static equations, or to otherwise highly specialized regimes aimed 
at very particular applications and model systems. In this paper we attempt to abate this absence by 
addressing the merits and limitations of a family of p-enrichment schemes implemented in a discontinuous 
Galerkin formulation to a generalized system of multicomponent equations. 

A brief note should be made here about our lexemic designation. That is, much of the literature on 
p-enrichment refer to it merely as p-adaptivity. This takes on a clear meaning, particularly in the context 
of /ip-adaptivity, where it generally describes what is happening in that context, in that p "adapts" along 
the lines of h in an often very complicated and highly hp-coupled manner. However, in the restricted 
case where one is only p-adapting, the context that emerges is a bit less subtle. Here, the general idea is 
to take some solution order p (where p = 1 would correspond to linears) , and attempt to raise the global 
average order of the solution without having to spend the computational resources required to raise the 
order of the solution in every element (e.g. p = 2 globally). Taken in this context it is then natural to 
view the ultimate goal of a purely p- adaptive scheme to be that of: global p-enrichment. Some authors 
however also prefer to use the term l p- refinement' instead of '^-enrichment.' Though this usage seems 
fairly clear in meaning, we prefer 'p-enrichment' to 'p-refinement' here, where the concept of refinement 
and coarsening seems naturally applicable to the mesh operations that occur in /i-adaptivity, while in 
the context of p-adaptivity we reserve enrichment and de-enrichment for the contrasting operations 
applied with respect to the order of the solution. Finally, it should be noted that P-enrichment may refer 
in chemical, biological and ecological settings to a type of phosphorous-enrichment, and should not be 
confused with p-enrichment in the numerical sense used in this paper, albeit in §4 to chemical problems 
involving the formation and depletion of aqueous phosphate levels. 

Within the context of /ip-adaptivity, the two volume work of Demkowicz, et. al. |15[ [16] provides 
both a nice introduction to some of the basic emergent features of Zip- adaptive regimes, as well as some 
advanced topics, such as /ip-adaptive schemes on boundary value problems applied to (a time-independent 
form of) Maxwell's equations; however, it should be noted that little emphasis is made specifically to 
discontinuous Galerkin (DG) formulations. Nevertheless, many of the basic principles that underlie 
/ip-adaptive methodologies can be viewed as relatively method-independent. For /ip-adaptive schemes 
applied specifically to DG formulations, the work of [45J applies an /ip-adaptive scheme to systems of 
nonlinear ODEs where an emphasis is made on exact numerical and mathematical behaviors of well- 
behaved solutions that demonstrate exponential rates of convergence in time. In [TJ] a nice /ip-adaptive 
scheme is presented which is based on using sharp a posteriori error estimates over (non) linear advective 
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PDEs. In [U [5] a parallel /ip-adaptive method is developed that uses a posteriori error tolerances in order 
to enforce sharp bounds on the error convergence. The /ip-adaptive scheme presented in [50] uses a pair 
of "smoothness indicators" over either the interior of the cell for p-enrichment, or with respect to the 
jump discontinuity between cells for /i-adaptation. In [10] an energy norm approach is given for Burger's 
and the Navier-Stokes equations, where artificial viscosity is added in the case of Burger's equations in 
order to develop the "diffusion scale" which serves as the implicit tolerance by which optimal convergence 
in norm is achieved. Additionally, we have recently developed an /ip-adaptive energy method in |38j 
based on a priori stability results, which use the local entropy of the fully coupled system of equations 
to "sense" and preserve global energy bounds, that we then apply to large multicomponent systems of 
reaction-diffusion equations. 

The results relating to p-enrichment, on the other hand, are more sparse. In the context of discon- 
tinuous Galerkin systems the two most prevalent results can be found in [9j 28J, which are discussed 
in detail in §3. Briefly, in [9] Burbeau and Sagaut develop what they refer to as "regularity sensors" 
for p-enrichment, and apply this model to Euler's equations in one and two dimensions from an applied 
engineering perspective. Next, the work of Kubatko, et al. in [28J expands on the model from [9], where 
a similar type of p-enrichment method is employed, but the I^-error behavior is carefully analyzed for 
dynamic p-enrichment with respect to varying fixed /i-refinement levels. In this work, both of these 
schemes would fit into the category of Type I fixed tolerance schemes, and both will be discussed in detail 
below in §3. 

We begin by introducing the generalized system of equations that we employ in §2, where we couch this 
system of initial-boundary value partial differential equations in the setting of a discontinuous Galerkin 
finite element method. Then in §3 we introduce and discuss a number of different general types of 
p-enrichment schemes, which include two types of both fixed tolerance as well as dioristic schemes. 

Finally, let us briefly discuss some of the engineering applications that we are interested in. Our 



generalized system of equations (see 2.1) encases a very large class of possible applications, but here we 
restrict ourselves to a pair of example applications arising in the field of coastal hydrology. In this setting, 
it has long been recognized that the shallow water equations offer a good and soluble approximation to flow 
regimes present in complicated physical systems (e.g. in coastal flows and storm surge models [6]. 113^ 1^9] . 
etc.). Here the depth-averaged shallow water equations (SWE) provide solutions to the elevation of the 
free surface of the water £ = x) and the velocity component u = u(t, x) of the flow as a fully coupled 
system of partial differential equations. Since this formulation models the essential mechanical properties 
of the flow pointwise as a function of time, it is clear that an advected scalar quantity (as we call t in §4) 
satisfying the linear transport equation 

Lf + U ■ V x l = 0, (1.1) 

coupled to the SWEs must track the flow characteristics of the solution, and thus for any chemically 
inert constituent t of the flow field represents the field of dispersal of that inert constituent, up to the 
depth-averaged approximation. Such a representation is particularly compelling when the constituent i 
is reasonably well-mixed (homogeneous) across the vertical stratification of the density field. Moreover, 
in such a case, it immediately follows that a chemical mass action principle must also be satisfied at the 
same corresponding scales. 

We present two separate models in §4 that employ this basic assumption in the context of coastal 
hydrology and the SWEs, and which attempt to take a first step towards addressing a pair of important 
questions concerning environmental coastal remediation. The form of the contaminant transport equation 
in §4.1 can be easily derived from first principles, though one should refer to [U fT2| \T§[ |4T>1 |4"T] for 
background on its uses and applications, and one should further note its close mathematical likeness 
to sediment transport models [39J. One thing to note is that we have suppressed the Fickian diffusion 



often present in the contaminant transport (1.1) aspects of the model. We have chosen this slight 
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simplification for this paper in part due to a frequent criticism leveled [30] against Eulerian frame solutions 
to the transport of passive tracers in contaminant models, which points out that numerical diffusion can 
become excessive for finite element solutions that are less than second order accurate in space. In fact, 
this provides one of the basic motivating factors for solutions that implement Lagrangian frame tracer 
particles in fixed lower order spatial schemes |31l l4"T] , which as a consequence demonstrate no diffusion 
at all (neither numerical nor Fickian), but rather travel indefinitely along the characteristic field. Since 
there remains a question as to what the role and relationship between the numeric dampening and the 
Fickian diffusion in the setting of a p-adaptive regime should be, where the polynomial order may span 
(by construction) the set p £ {0, 1, . . . , n}, we save this topic of macroscopic particle diffusion for another 
time (though we also refer the reader to |38| for some preliminary results in this direction). 

Our second example in §4.2 takes the standard transport model from §4.1 and extends it to include 
reactive chemical constituents. We do this by coupling a standard chemical mass action term. This type 
of reactive shallow water model has been implemented recently in [21 [11] , and can be inferred from many 
related systems, such as [U [32]. We provide a fully Eulerian solution to the multicomponent reactive 
system of equations, by implementing the discontinuous Galerkin model for chemical reactor systems 
developed in |38] , and apply this system to a realistic mesh of the Brazos estuary emptying into the Gulf 
of Mexico, and analyzed in the context of the development of contaminant induced "dead zones" in these 
types of regions. 



§2 Governing equations and formulation 

We are primarily focused in this paper on systems related to the shallow water equations [6j , though in this 
setting we wish to include systems that have both complicated reaction dynamics and free boundary data. 
Thus, we generalize somewhat to consider the following initial-boundary value problem in Q x (0, T), 
taking f2 C M 2 with boundary d^l, such that the system is governed by: 

Ut + F x — G x = g, given initial conditions U\ t=0 = Uo, (2-1) 

with Robin free boundary, 

mUi + V x U i>x (h ■ n + q • r) - = 0, on dtt - (2.2) 

More precisely, the system is comprised of an m-dimensional state vector U = U(t, x) = (Ui, . . . , U m ), 
an advective flux matrix F = F(U), a viscous flux matrix G = G(U ,U X ), and a source term g = 
g(t,x) = (gi, . . . ,g m ), where x £ R 2 and t £ (0, T). The vectors a, 6, c and / are comprised of the m 
functions, en = ai(t,x), hi = bi(t,x), c\ = Ci(t,x) and fi = fi(t,x) for i = 1, . . . , m, where n denotes 
the unit outward pointing normal, and t the unit tangent vector. The time-varying free boundary 
dQ, f ree = d£lf ree (t,x) adjoined to the absolute domain boundary dQ make up what we refer to as "the 
effective boundary," which we denote by d£lo = d£l U d£lf ree (for example, we use a wetting and drying 
boundary condition below that corresponds to the free boundary). 

In addition, because we are interested in approximate numerical solutions of the form of [HE] restricted 
in part to the family of LDG (local discontinuous Galerkin) methods as developed originally for elliptic 
equations [3J, we rewrite ( |2.1[ ) as a coupled system in terms of an auxiliary variable S, such that the 
system we solve becomes: 

U t + F x -G x = g, and S = U x , (2.3) 

where we have substituted the auxiliary variable into the viscous flux matrix, so that G = G(U, S). 

For notational completeness we adopt the following discretization scheme motivated by [171 [37]. Con- 
sider the open set !!cR 2 with boundary dil, given T > such that Qt = ((0, T) x Q). Let ^ denote the 
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partition of the closure of the polygonal triangulation of O, which we denote f^, into a finite number of 
polygonal elements denoted Q e , such that = {O ei , £l e2 , . . . , 0, ene }, for ne G N the number of elements 
in Qh. In this work we define the mesh diameter h to satisfy h = min.ij(dij) for the distance function 
) and elementwise face vertices Xi,xj G d£l e when the mesh is structured and regular. For 
unstructured meshes we mean the average value of h over the mesh. 

Now, let Tij denote the face shared by two neighboring elements Q ei and fi e ., and for i £ I C Z + = 
{1, 2, . . .} define the indexing set r{i) = {j G I : O e . is a neighbor of £l ei }- Let us denote all Q ei containing 
the boundary dilh by Sj and letting Is C TLr = { — 1, —2, . . .} define s(i) = {j € Ib ■ Sj is a face of Q ei } 
such that = Sj for S7 6i £ fi^ when Sj G <9f2 ei , j G Ib. Then for Sj = r(i) U s(i), we have 

do e ,= (J Ty, and 9n ei naf2 h = |J rv 

j'eE(i) ies(j) 

We are interested in obtaining an approximate solution to [/ at time i on the finite dimensional space 
of discontinuous piecewise polynomial functions over restricted to 5^, given as 

s p h (n h , sr h ) = {v : g &> p {n ei ) vn ei g sr h } 

for ^ p (f& ei ) the space of degree < p polynomials over Q ei . 

Choosing a set of degree p polynomial basis functions N( G <^ p (f2 ei ) for £ = 1, . . . , n p the corresponding 
degrees of freedom, we can denote the state vector at time t over £l ei , by 

Tip 

E7 fcp (t,a!) = 5^l7j(t)iyj(x), Vxefi e „ (2.4) 
e=i 

where the NVs are the finite element shape functions in the DG setting, and the U\ J s correspond to the 
unknowns. We characterize the finite dimensional test functions 

Tip Tip 

v h p,Wh P G W Kq (9, h , & h ), by v hp (x) = J^t>jiyj(x) and w hp (x) = Y J w e N e( x ) 

1=1 £=1 

where v\ and w\ are the coordinates in each fi ei , with the broken Sobolev space over the partition ^ 
defined by 

W k >i(n h , <%) = {u : wm.. G ^ fc '«(Q e J Vft e , G 3T h }. 



Thus, for U a classical solution to (2.3), multiplying by v^ p or w^ p and integrating elementwise by 
parts yields the coupled system: 

— / U v hp dx + I (F ■ v hp ) x dx - [ F : v hp dx 
dt Jn e - Jn e . Jn e . 

°z L i L j 

- / (G ■ v hp ) x dx + G : v^dx = v hp - gdx, (2.5) 
Jn e . Jn e . Jn e . 

c ! c ! c i 

/ 5] • Whpdx — I (U ■ Whp) x dx + U : w^ p dx = 0, 

where (:) denotes the scalar product. 

Now, let rijj be the unit outward normal to dQ ei on T^, and let «|p.. and U|r- 4 denote the values of v 
on T^ considered from the interior and the exterior of Q ei , respectively. Then by choosing componentwise 
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approximations in (2.5) by substituting in (2.4), we arrive with the approximate form of the first term of 

(2.6) 



(2.5) given by, 



d 
dl 



I U hp - v hp dx « f U v hp dx, 
the second term using an inviscid numerical flux 3>j, by 

^i{U hp \T 13 ,U hp \ Tji ,v hp ) = / ^{UhpWi^U^r^riij) -Vhplr^dl 

,-6B(<) Jr *i 



f 2 

jes(i) i=i 



(2.7) 



and the third term in (2.5) by, 

@i(Uh P ,v hp ) 



F hp : v!*da 



F : v h fdx. 



(2.8 



Next we approximate the boundary viscous term of (2.5) using a generalized viscous flux & such that, 

{^h P ,u hp ,v hp ) = y i 

(2.9) 



while the second viscous term is approximated by: 



<A{T. hp , U hp ,v hp ) = G hp : v hp dx « / G : v hp dx. 

For the auxiliary equation in (2.5) we expand it such that the approximate solution satisfies: 
£!i(U,Yl hp ,U h p,w h p,Wx P ) = / T, hp ■ w hp dx + / U hp :w^. p dx 

^2 / 17(17^1^., l/fcplr^jtfAplry.nij)^ 



(2.10) 



(2.11) 



where, 



YY U(U h p\ Vi .,U h p\ Vji ,w h p\r iV n ij )d'E » ^ / ^(E/)z • (n^M^li^d: 
ie/ iesfi) r ^ is/ jes(i) 1=1* 



given a generalized numerical flux [/, and where 



Sftp • Wh p dx w / I] • Whndx, and / [//,,„ • w^Pd. 



'hp 



hp 



Jn, 



U ■ w* p d 



x. 



Combining the above approximations and setting = ^Q e . g ^ while defining the inner product 

i a h P i bh P )n g = / a hp " 



we arrive at our approximate solution to (2.3) as the pair of functions (U hpi ^hp) f° r an * £ (0>^) 
satisfying: 
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The Discontinuous Galerkin formulation 

a) U hp GC\[0,Ty,S p h ), S ftp e^ 

b) ^ (U hp , v hp )a g + ®(U hp , v hp ) - @(U hp , v hp ) 

- &(E hp , U hp , v hp ) + J^CE hp , U hp , v hp ) = 0, ( 2 - 12 ) 

c) .2(17, E hpj U hp , w hp , w%>) = 0, 

d) U hp {0) = U hp U , 

where Uh p is a projection operator onto the space of discontinuous piecewise polynomials S?. Below 
we utilize the standard L 2 -projection, given for a function / G L 2 (J7 e .) such that our approximate 
projection f Q h G 7 2 (^ e J is obtained by solving, J Q f oh v hp dx = j n f v hp dx. 



The discretization in time follows now directly from (2.12), where we employ a family of SSP (strong 
stability preserving, or often "total variation diminishing (TVD)") Runge-Kutta schemes as discussed in 
|444 146j . That is, for the generalized SSP Runge-Kutta scheme we rewrite ( |2.12^ ) in the form: ME/t = R, 
where U = (U\, . . . ,U P ) for each element from (2.4), where R = R(I7, E) is the advection-diffusion 



contribution along with the source term, and where M is the usual mass matrix. Then the generalized s 
stage of order 7 SSP Runge-Kutta method (denoted SSP(s, 7) or RKSSP(s, 7)) may be written to satisfy: 

17(0) = jrjn 

i-l 

E/W = {<XirU r + AtfiirM^BT) , for i = l,...,s (2.13) 

jjn+l = c/ ( S ) j 

where R r = R(C/ r ,5] r ') = R (U r , S r , t n + 5 r At) and the solution at the n-th timestep is given as 
U n = U\ t=t n and at the n-th plus first timestep by U n+1 = C7i t=t n+i, with t n+l = t n + At. The a. LT and 
(3i r are the coefficients arising from the Butcher Tableau, and the third argument in R r corresponds to 
the time-lag complication arising in the constraints of the TVD formalism. That is 8 r = ^4=0 f^rl 1 where 
Hir = Pir + Sz=r+i Plr&ih where we have taken that a{ r > satisfying ^r=o air = ^ , 

Our examples in this paper will all be given in the context of the discontinuous Galerkin shallow water 



code described in [HI [131 125H27] [29] . which employs a fully coupled system of (2.12) including coupled 
eddy viscosity, time varying free boundary data, coupled chemical reactor models, etc. to the shallow 
water system of equations. For the polynomial basis we choose the hierarchical Dubiner basis, and our 
meshes are comprised of triangular elements. 



§3 Types of ^-enrichment 

Here we present a family of generalizable dynamic p-enrichment schemes based on both local and global 
data. Our first class of methods effectively extend the formalisms presented in [36J to what we refer 
to here as fixed tolerance schemes. These schemes use only local data to estimate the regularity of the 
solution, and then adapts the solution based on a fixed scalar-valued global tolerance setting. Next we 
introduce the class of dioristic schemes. In these schemes, the solution is adapted based on global bounds 
on the variation which are chosen with respect to properties of the local regularity. Here, the smoothness 
estimators from |36j are extended to treat the solution vector Uh p as a conjunction of weakly coupled 
components which each have their own variational bounds with respect to a set global tolerance. This 
variation in smoothness is then used to determine the stabilization regime. 
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3.1 The fixed tolerance approach 

Let us first present the class of dynamic-in-p enrichment schemes that rely on the local solution data 
and are characterized by a global fixed tolerance. Here we may state our general goal as being one of 
two things: (1) to locate within the solution areas of "excess" variation and to probe these areas in order 
to capture higher order structure, or (2) to identify areas of potentially destabilizing variation in the 
solution, and to delimit these areas for de-enrichment in keeping with the prescribed stability conditions 
of the solution (e.g. the CFL condition). We will discuss these more below. 

The first type of enrichment scheme, which we refer to here simply as the Type I fixed tolerance 
scheme, apply to solutions in which substantial regularity (e.g. U G C°°) might be assumed over the 
entire domain f2 x [0, T]. In these situations, the regularity of the solution is to be exploited such that 
we aim to resolve the "areas of highest interest," which are those areas which show maximal variation of 
the solution U hp- The Type I scheme is particularly attractive in application models where local high 
energy behavior is the signature of an "area of interest." For example, in coastal models of hurricane 
storm surge, where one might want to resolve wave models along bay and river inlets in order to recover 
high accurate flooding behavior, this type of scheme might be particularly well-suited, as narrow channels 
often tend to focus the relative energy signature of a local subregion. 

To be more precise, the Type I scheme takes the approximate solution vector U hp and computes an 
auxiliary sensor II (which we may also think of as a "relative smoothness" or "relative regularity" sensor) 
over each i-th component of U hp (where i = 1, . . . , m), defined by: 



IT 

j 



TP I 

u hp\ 



TP I 

u h P \ 



Xj 



(3.1) 



where the solution Uh P is evaluated at c, the centroid of element f2 e , and LVj, the midpoint of the j— th 
face of f2 e . The function Xj m ay be set to either the distance Xj = ~ c l as i n [28| , or the product 
Xj = ojjC as in [9]. In either case, over each timestep n with respect to (2.13) the following p-enrichment 



functional (£i e = <£i e (£P (f2™)) is evaluated over each cell f] e 





Type I fixed tolerance scheme 






' <? k+1 (n n e ) if ((sup, sup, H] > e) A (k + 1 < p max )) A (r > t w ), 






@> h -\V%) if (inf, sup, 11} < e) A (k - 1 > p min ) A (r > t w ), 


(3.2) 




^(£1™) otherwise, 





where To is a counter that restricts the enriching/de-enriching so that it only occurs every t w £ 
timesteps, and where k E {1, ... ,p}. 

Depending on the choice of the global tolerance e — which is just some positive number e S M + 



N 



the Type I fixed tolerance scheme (3.2) can provide either a fairly stringent or a fairly loose cutoff with 
respect to which elements it flags for p-enrichment. The most immediate difficulty that the scheme seems 
to present, is how exactly to choose the value of e. If too small, then the entire solution immediately 
climbs to p m ax and the dynamic nature of the algorithm is lost on the domain, and nothing particularly 
useful is accomplished. If the value is too large, then the solution gets trapped at p m i n and never 
enriches appreciably beyond the initial p state of the solution. In fact, this issue becomes somewhat 
of a complication for solutions which demonstrate substantial variation over time with respect to their 
"regularity profile." In this case, the question arises: is it better to tune the value of e to the initial state 
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of the solution?, or is it better to tune e to the maximal state of the solution?, or is it better to tune e to 
the average state of the solution over (0,T)?, ... and so forth. 

It should not come as a surprise here that the answers to many of these questions are quite application 
dependent. However, the questions themselves reveal perhaps the most important motivation for looking 
into p-enrichment regimes which do not harbor this type of global restriction with respect to the fixed 
tolerance scalar e. Admittedly, one could also work to generate an application dependent form of e which 
varies appropriately with respect to the solution, as e = e(t,x), in order to stay "centered" with respect 
to the average value of IT*-, for example. And while such an approach may actually be quite effective for 
a specific model, it would also be quite difficult to generalize without addressing in some detail the form 
of the model-dependent regularity estimator itself. We will revisit this issue some in §3.2. 

The second type of enrichment scheme we would like to consider in this section is the Type II fixed 
tolerance scheme. Here, in the Type II regime, we omit any assumptions on the regularity of our solution 
a priori. Rather, in the Type II regime we assume the possibility of numerical shock fronts arising due 
to instabilities that might be caused by pathological perturbations in highly variable solutions, natural 
discontinuities developing in the nonlinear systems, artifacts arising from numerical instabilities, or the 
like. 

One way of formalizing when it is appropriate to choose such a regime, is to say that Type II schemes 
are reserved for solutions demonstrating both appreciable local gradients V x Uhp ^ 0, and nonzero local 
mean curvature of the solution, related by ~ V x ■ V x Uh P (see [13)- As such they are designed with an 
eye towards solutions which have both a large dynamic range in the magnitude of the unknowns, as well 
as solutions with large local spatial variations. 

In order to fully develop the Type II schemes, we again develop a regularity functional for Uh P , but in 
this context we employ a slightly different formalism from the Type I schemes. That is, the Type I fixed 
tolerance schemes abstractly try to capture the amount a solution varies with respect to the element's 
center and the members of its boundary. But now, we wish to try and develop an auxiliary functional 
n*- which takes into account the order of the variation over the element. 

Perhaps the most obvious way of developing such a functional would be to utilize the jump condition 
between the base element and its neighbors. After all, if the solution varies dramatically between two 
cells in a discontinuous basis with respect to the jump condition, then clearly the order of the variation 
is high in that area. For example, in [9] the Van Leer minmod function is utilized to define the auxiliary 
sensor 

Hj = minmod(E/y v + - U% p \ c , U\ p \ v - - U\ p \ c ), (3.3) 

where Vj is the j-th vertex of Qh p evaluated with respect to the standard jump condition. As — > 
the solution becomes smoother, and one may subsequently employ an algorithm similar to ( |3.2[ ) (up to 
the direction of inequalities, for example). 

This method offers a nice solution to the problem, but it introduces two factors which we wish to 
avoid in our present context. The first is, it introduces a nonlocal stencil so that the auxiliary functional 



(3.3) does not strictly depend on only the information contained within a local element. Which leads 



to the second issue, being: the auxiliary functional (3.3) would then be most effective when the mesh 
itself is aligned along the lines of maximal variation. But this is precisely what occurs in many standard 
/i-refinement algorithms (e.g. see |15|. I38j). which use the solution behavior across element boundaries in 
order to refine the mesh along potential discontinuities — a feature considered a requirement in order to 
obtain the requisite exponential convergence of hp- adaptive regimes. 

However, in the context of an /ip-adaptive regime trying to /i-refine and p-enrich a cell simultaneously 
can rapidly lead to unstable behavior. This is impossible to avoid when using the same regularity indicator 
for both h and p, unless one adopts one of two alternative approaches: either /i-refine the cell while p- 
de-enriching it — leading to an automatic counteraction in the accuracy of the local solution, which is 
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unnecessary — or, choose two different tolerances e for the h- refinement and p-enrichment, respectively. 
While this may work in practice, it erroneously seems to imply that the role of h and p are exactly the 
same, but just "act" at different magnitudes. Since this is not the case (particularly at discontinuities), 
we look to develop a measure which isolates the unique properties of p-enrichment more explicitly. 

Towards this goal, we define a local (cell-dependent) regularity estimator that relies on the weighted 
norm of the higher order components of the solution. That is, let Uhp be the evaluated solution in the 
basis where the highest components have been truncated. In other words the members of & >k ~ 1 (fl2) in 
our hierarchical basis. Then we take the norm of the higher order components (i.e. all those coefficients 
corresponding to A; > k— 1), and quotient out by the norm of the solution at k, defining the local regularity 
estimator: 

/ _ II \ 

n f= trri „ MM(ne) , for U\ p G ^> k ~ l {W e ) and U\ p E & k {W e ). (3.4) 

y \\ U hp\\L1{p.e) J 

Here we have taken the standard Lebesgue L q norms (except when q = 2 in which case we preferentially 
take the standard inner product). 

Now our regularity estimator ( |3.4[ ) can be viewed as a weighted average of the nonlinear coefficients 
of the basis (at least whenever p m in = !)• We use this function then to develop our Type II algorithm, 
which is designed to truncate higher order oscillations by way of p-de-enrichment, while p-enriching areas 



which show high regularity. That is, in the Type II fixed tolerance scheme, we replace (3.2) with the Type 
II functional: 





Type II fixed tolerance scheme 






if (su Pi log 10 n? < A) A (k + 1 < p max ) A (to > t w ), 




<£n e = { ^ fc - 1 (^; 


if (infi log 10 n? > A) A {k - 1 > p min ), 


(3.5) 




otherwise, 





where the bound satisfies 

A = / logl o 5k ~ 9 + c ' for P > P mi " 



supjlog 10 nf, otherwise 



such that c, c S M + are user defined constants, where c E (0, 10) is recommended (see for example |50j ) 
for resolving discontinuities in the context of /ip-adaptivity, and where we have found c 6 (—2, 2) optimal. 



The basic intuition that underpins the use of (3.4) is the observation that the coefficients in the basis 



are assumed to decay at a rate comparable to that of the Fourier coefficients in a standard expansion of 
the solution — which clearly decay at a rate of 1/fe 4 for q = 2 (see |42| |4~3"1 150j). to which we obtain an 
indicator of the relative local regularity of the solution, i.e. the faster the coefficients decay, the more 



regular the local solution. Thus we obtain equation (3.4), which approaches zero as the solution becomes 



smoother, and where setting c > is a sharper restriction than the more permissive (though substantially 
less stable) c < 0. 

The Type II fixed tolerance scheme can be very stabilizing due to the fact that given an appropriate 
choice of A it will always de-enrich with respect to the elements experiencing the maximal variation. 



However, as with the Type I fixed tolerance scheme (3.2), finding the correct value for A can become 
a very subtle procedure, and the choice may ultimately be far from ideal over the full solution domain 
0, x [0, T]. Moreover, the Type II regime may not satisfy the objective of the enrichment scheme in the 
context of the given application, where one might be more concerned with fleshing out structure in areas 
of higher variation, rather than optimizing the algorithm to enrich most stably. 



3 Types of p-enrichment 
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Nevertheless, the existence of solutions which may vary substantially over either space or time (or 
both) underscores the need to develop schemes which do not rely on globally fixed space and time 
independent constants. Though a number of approaches have been developed along these line, here we 
present a novel set of algorithms which inherently rely upon a "stabilizing center" condition, and to which 
we refer to conventionally as: dioristic enrichment schemes. 



3.2 The dioristic algorithms 

The class of dioristic enrichment schemes do not rely on knowing any global properties of the solution 
ahead of time (e.g. knowing the global bounds on the local variation of the maximal and minimal 
components in order to choose an appropriate tolerance such as e or A, as required in §3.1), but is designed 
to maximize the available computational resources in such a way as to distribute as many enriched cells 
as possible with respect to a parsing of the total local variation present within the components of the 
state vector around a "stabilizing center." In order to achieve this, we must define what is meant by 
a "stabilizing center," as well as develop a functional representation of our solution which couples each 
component of the state variable in a way that consistently identifies the behavior of the local variation. 

A particularly nice way of developing such a fully coupled functional over Uh p is to derive a unique 
energy functional 5? = -5^(t, x) for the system. However, deriving such a functional 5? requires both 
choosing the exact form of the system a "priori (often including, for example, fully explicit forms for 
the boundary conditions) as well as knowing quite a lot about the form of that equation (e.g. analytic 
existence and regularity of solutions). Thus, these types of methods are far more difficult to generalize. 
As an alternative we wish to work in a more general framework, where we assume that the formal entropy 
of the system may not be easily available or implementable a priori, and so we resort to a slightly more 
weakly coupled form of the energy functionals. Let us refer the reader to [38] for more details on how to 
construct entropy(,5^)-stabilized /ip-adaptive regimes by way of formal energy methods. 

We now proceed by defining what is meant in the present context by a "stabilizing center." Since our 
goal is to develop an efficient and easily generalizable p-enrichment methodology, we assume from here 
forward that we may only use the information at timestep t n from either U^p or the derived regularity 
estimators H(Uh P ) in order to determine our dioristic functionals. To do this we first define the range 
£ = £(n.(Uhp)) of the variation in the regularity of each component i of the solution vector Uh p as 
£i = (maxQg Hi — minQ g IT). That is, by the global maximum maxQ g (-) we simply mean max^g = 
maxvn e en g (•)> an d so forth. Next we set the composite function Si = Si(£i, fa) of the range £j and the 
adjustable weight fa £ (0, 1), such that the product S{ = fa£i determines a weighted distance of the range. 

Then denoting the global average smoothness of in each component i of the solution vector as Avg^IT, 
we define the "stabilizing center" at timestep n as the discrete subdomain c C Q,^ p comprised of the union 
of elements over which the solution satisfies the condition, 




IH- Avg^IIil < Si, Vij. (3.7) 
The subdomain c is "stable" in the sense that the mean value of the regularity of the solution must 



be commensurate to the stability settings of the free parameters associated to the formulation of (2.12) 



and ( |2.13 ). Informally what this means is that "on average" the free parameter settings (e.g. At, dx, 



etc.) must be chosen such that the stability conditions (e.g. the CFL condition, etc.) are "on average" 
satisfied. In fact, this should be viewed as an implicit condition on both of the dioristic algorithms, which 
can be stately more heuristically as simply assuming that the stabilizing center c is in fact stab le. 



Then we are able to determine the Type I dioristic functional as complementary to (3.2). That is 



here, as in §3.1, the Type I functional aims to enrich areas of greater variability using IT from 



3.1). But 
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in the dioristic setting, the enriching about a fixed global parameter e is replaced by an enriching about 
the "stabilizing center" c. That is, over each component i of the state vector Uh p the Type I dioristic 
functional 2)^ is set to satisfy: 





Type I dioristic scheme 




' <^ k+ \n2) if 


(| su Pj U) - Avg Qg su Pi nj| > S t ) A (fe + 1 < Pmax )) A {t > t w ) , 


= < 


&> k - 1 {^l) if ( 


| su Pi u) - Av gng su Pi n* | < St Vi) A (fe - 1 > p min ) A (r > n, 






otherwise, 




(3.8) 



where we are using the same definitions as in §3.1, except that here Si has the following functional depen- 
dencies: Si = 5i{supj IPj, Hi). More precisely, in this context we simply define Si = ^(max^g, sup^ II* — 
minn g sup,,- II*) for each component % of the state vector, where the global average smoothness takes the 
form, Avgj^ sup^n}. 

The Type I dioristic scheme may be described as an algorithm that takes the value of each component 
of the solution at timestep n, finds how close to the average value of the component over the whole domain 
the evaluated solution is; and then, when any component of the solution's variability exceeds that of c 
on the element, it locally p-enriches the solution. On the other hand, the Type I dioristic functional 
only de-enriches the solution if the element lies within c, thus having the effect of reducing the average 
polynomial degree in c. This feature is developed to attempt to counterbalance the destabilizing effect 
that p-enriching the solution only in the areas of highest variability can have over even a relatively uniform 
space of solutions. 

For the Type II dioristic functional we develop a similar approach contextualized with respect to 
the "stabilizing center" c. The aim of the Type II dioristic scheme is to develop a scheme which takes the 
elements within the stabilizing center c and elevates their polynomial order, while flagging those elements 
in the extremal regions {i.e. Vt ei ^ c) of for de-enrichment. 

Here again we use the composite function Si for each component i of the state vector Uhp, but now we 
define the global average regularity of each component of the solution by, Avg^ g log 10 II?. Then similar 
to the Type I dioristic functional, we define: 



Type II dioristic scheme 

& k+1 {K) if (| logio nf - Avg ng log 10 n?| < S t Vi) A {k + 1 < p max ) A (r > t w ), 
2>n e = <j & k - l {W e ) if (| log 10 nf - Avg ng log 10 n?| > S,) A {k - 1 > Pmin ), 

£? k {Q r e ) otherwise, 

(3.9) 



where now we have functional dependencies given by Si = Si{Uf, Hi), and defined by Si = /ij(max^ e log 10 II?- 
mm Qg log 10 II?). 

The Type II algorithm obeys a certain type of parsimony with respect to its enriching functionality, 
where the only way for an element to get enriched is if the solution vector is smooth enough on the element 
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to become a member of stabilizing center c. Then, within these islands of stability the polynomial order 
is increased. Everywhere else in the domain, the solution is de-enriched, operating under the tacit 
assumption that elements Q ei ^ c possess either potentially destabilizing irregularity, or, relative to the 
rest of the solution, are "nearly constant" so that the higher order structural information recovered 
represents a meager benefit to the global behavior of the solution. 

Notice that both the Type I and Type II dioristic algorithms require setting the local parameter fa. 
Since the parameter 5i in both (3.8) and (3.9) is determined with respect to the range the parameter 
fa just represents the weighted distance with respect to the total variation of the regularity at timestep 
n that the cell based regularity may vary from the global average regularity in order to remain a member 
of the stabilizing center. It is also possible to allow fa to depend on time fa = fa[t), though due to the 
strong time-dependence of £j = £i(t) this is usually unnecessary in practice, and, as the parameters in 
§3.1, the behavior of fa in time is quite subtle to try to accurately predict. 

Nevertheless, both dioristic algorithms use the componentwise local variation in the solution space 
with respect to the the global variation in that component at a particular time t n to determine which 
areas in the domain are — relatively speaking — experiencing the largest relative fluctuations. Thus for 
well-behaved solutions the dioristic functionals £>j ! and can be interpreted as energy-type functionals 
which sense the local variation in each component of the solution, and then p-enrich the domain with 
respect to some weighted percentage of the relative variation in each component as determined by a 
choice of fa, while trying to stabilize by simultaneously de-enriching other components of the solution 
with respect to a "stabilizing center." 



§4 Some numerical examples 

We now consider a number of examples in order to test and analyze the behaviors of the various p- 
enrichment schemes presented in §3. The first example is designed as a difficult problem with many 
symmetries, where the variation of the solution is quite large, existing at the precipice of numerical shock 
formation. The second example is quite complicated and is more of an application model, which attempts 
to realistically model estuary eutrophication in the Gulf of Mexico. 



4.1 Contaminant plume declinature 

Our first example is that of a plume declinature model. The underlying physics that the model tries 
to address — with respect to the fairly idealized setting in two dimensions — is, consider a bay with 
a contaminant plume on the far side of an inlet into the bay (see Figure [I]). Then we ask the physical 
question: is it possible to prevent the contaminant plume from entering the bay by imparting mechanical 
wave energy on the near side of the the bay using a standard hydraulic wavepool generator? As we 
show below in two dimensions and given a very powerful local hydraulic generator at least, the answer 
is yes. We additionally address the issue of p accuracy in the p-adapting regime, versus the standard 
convergence-in-p of the solution in a slightly more restrained setting. 

For the solution to the contaminant declinature problem, consider the transport form of the two- 
dimensional shallow water initial-boundary problem, determined by: 

d t ((i) + V x (Hlu) = 0, 
d t (Hu) + V x e + S- v^x(Hu) - g(V 

6 = (hu ® u + \g{H 2 - h 2 ) 



r h = 0, 



(4.1) 
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Figure 1: Initial conditions, boundary conditions and course mesh of the contaminant plume declinature 
model, where the bay is on the right with a wavepool forcing. 

where the initial data satisfies 



t=0 



Co, Ut=Q = Uq, L\t=0 = t0- 



Here the solution space is comprised of the velocity field u = u(t, x), the elevation of the free surface 
( = x), the total water column height H = Q-\-b, where b = b(x) is the time independent bathymetric 
depth (measured positive downwards), and i = i(t, x) represents the relative concentration of a chemically 
inert contaminant. The terms containing the gravitational constant g E M + correspond to those deriving 
from the frequently employed hydrostatic pressure assumption, while the constant rj £ M + is the eddy 
viscosity coefficient (here set to r\ = 10 m 2 -s -1 ), and S corresponds to the remaining source terms, which 
generally could contain Coriolis forces, bottom friction, wind forcing, etc. Here we simply set a bottom 
friction approximation using the Chezy formula such that S = .0025|tt|/i? (see §4.2 for more complicated 
source settings). 

The initial data is Co = m, b = 10 m, and uq = m-s with initial relative contaminant concentra- 
tion: 



to 



1/5, for < r < Oi, where r = yjx + a) 2 + 



420 m, and a\ = 250 m. 



It is not difficult to see that the system (4.1 ) is easily formulated in the form of (2.3), where here we use a 
local Lax-Friedrichs flux, and where we use an upwinding scheme for the coupled contaminant transport 
model. 

The boundary conditions are separated out so that we have a union of Dirichlet conditions over the 
total domain boundary: dQh = d^sand U d^tide U deforce- Each component of U^p is set independently, 
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Table 1: We give the normalized L 1 -error of ih p after half a day {i.e. T = 0.5 days) of simulation time with 
respect to the p = 7 converging in p solution using af orce = m, with the relative CPU times (rCPU) of 
the approximate solutions over continuously adapted (e.g. p = 1-3 implies that p spans {1, 2, 3}) solutions 
coupled to the the enrichment schemes. 
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Figure 2: Here we show the p-convergence of the solution denoted in Table [TJ where the static forcing is 
used. 

but due to the physical constraints on the system can be simplified to satisfy, 

Cb = (Cb,l)sand U {(betide U ((b,3J force, ^b = ( ^6, 123 ) sand, tide, j or ce , 

Ub = («6,n,l • n + u b,r,l ■ T)sand U (u&, n ,2 " " + W6,t,2 " 1~)t«fe (4.2) 

orce ■ 

The first boundary type dVL san d corresponds to a sand beach boundary condition, the second dfltide to 
an open ocean tidal inlet condition, and the third d£l force to either a hydraulic wavepool periodic forcing 
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Table 2: We show the robustness of the p-enrichment, as well as the effects on the accuracy of the weak 
entropy in strong boundary forcings. Here we show the normalized Z^-error of ih p after 2 days (i.e. T = 2 
days) of simulation time with respect to the p = 7 solution using a/ orce = 6 m, with the relative CPU 
times (rCPU) of the approximate solutions. Here * means the solution is run using an adapting-in-p 
slopelimiter (discussed below). 



condition, or a static condition (see below for details). Here, as in (2.2), r is the tangent unit vector and 
n is the outward pointing unit vector. 

Spatially, the entire west edge of the domain as labeled in Figure [T] is comprised of a tidal inlet 
condition delude along 20 elements, and is 3000 m long. The center of the east edge of the bay — which 
is defined as the two east edge elements with \y\ < 150 m — is given the forcing condition d£lf orce . All 
of the remaining 92 edges of the domain boundary dVth are given the sand beach boundary condition 

Let us define these three different boundary types precisely. The first order transmissive scalar 



boundary values from (4.2) are determined by the value on the elements interior at the boundary, so that 
(b,i\dn e . = (\dtt e . an d tb.2\8Qe. = L \dn e - The scalar forcings at the boundary, on the other hand, are set 
so that for I G {1,3} they satisfy: 



(b,2\dn ei = a t id e cos(uJtidet), (b,3 \ dQ ei = "-force \ cos(w/ orce t) | , and ib,l\dn H =Q- (4.3) 

The amplitude and frequency of the tide are given as aude = 0.75 m and uitide = 7.02 x 10~ 5 rad • s _1 
respectively, while the amplitude and frequency of the wavepool forcing is given by cif orce = 6 m and 
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w force = 2.56 x 10 -3 rad-s -1 respectively, while the static forcing is determined uniquely by a/ orce = m. 
Both amplitudes are ramped up over half a day in order to avoid shocking the system, where a tide of 
three quarters of a meter occurs every half day over dVttide, an d a very large and localized wave of six 
meters is generated every twenty minutes over dQ.f orce with respect to the bay (see Figure [TJ. The third 
condition on the relative contaminant concentration in (4.3) represents first a filtration condition at the 
wavepool inlet, and second an absorption assumption at beach boundaries {e.g. petroleum contaminated 
sand HU). 

Finally, the velocity boundary conditions are set to obey a dampening with respect to the sandy 
beach and an open ocean quenching condition, while the wavepool vector is chosen as to quench the 
nonlinear velocities forced by Q\dn H - More clearly, for I G {1,2,3} we have that Ub !n ,i\dn ei = m • s _1 
and u b>Tj i\ d Q e , = m • s _1 . 

Let us briefly discuss the results of our numerical experiments. First, it should be clear that our 
boundary forcings drive our solution, and are weakly imposed as discussed above. These weak entropy 
conditions have been analyzed in previous work |34| [3"7] , and are known to be a source of error in the 
approximate solution. In fact, the error on the boundary elements may scale with the mesh size h in 
our adapting regime, such that in our case, where are boundary forcings are very strong relative to a 
small interior domain, the weak entropy boundary error frequently dominates. Because of this, we do 
not achieve (or expect to achieve) the rates of convergence in p one might expect from a manufactured 
solution, etc. in the case of the wavepool forcing, which is quite strong. In fact, we push the boundary 
forcings to the edge of what the convergence- m.-p can even distinguish, as is clear in Table [2] between the 
p = 7, p = 5 and p = 1 cases; where our focus is rather the robustness of the enrichment schemes. In 
Table [2] we see impressive robustness of the Type II algorithms, even with rigid constraints imposed by 
both the CFL condition and on the weak entropy boundary forcings, where it is clear that even running 
the solution in quite dynamic circumstances (e.g. when p £ {1, . . . , 7} and the wavepool forcing on the 
boundary is applied) the solution is still stable at quite large timesteps (e..g At = 2 seconds) without 
introducing any observable loss in accuracy. 

In contrast, when we suppress the strong boundary layers given by the wavepool forcing, and employ 
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Table 3: Here we provide the p value as computed over each element corresponding to the four p- 
enrichment algorithms from §3 given the hydraulic wavepool forcing. Again * corresponds to the use of 
the slopelimiter from |36j 

the static forcing condition a/ orce = m instead, the weak entropy on the boundary is reduced and 
propagates below the levels of the accuracy of the interior solution. We present these results in Table [TJ 
Here we see the requisite p-convergence behavior of the solution with the reference solution at p = 7 (as 
shown in Figure [2]). Notice that all four schemes under the correct settings may be used to approximate 
p-accuracy lying between integral values. From the timestepping chosen, it is immediately clear that 
the Type II algorithms demonstrate greater relative stability compared to the Type I algorithms. In 
fact, it should be noted that the Type I algorithms require settings that are either largely weighted 
towards p = p m i n or p = p max - In other words, in this test case, because the Type I algorithms enrich 
in the areas of greatest variation, in order to achieve stability the amount of p adaptation over T must 
be minimized. In contrast, the Type II algorithms demonstrate substantial robustness at a number of 
different timesteps and settings, while still achieving p-convergence. Nevertheless, it must be stressed that 
in some applications one strongly desires structural information in areas of maximal variation; so much 
so that the loss of stability in the enrichment process pays for itself in terms of recovering information 
with respect to the proper locale of interest. 

The different behavior of the fixed tolerance algorithms compared to the dioristic algorithms is also 
emphasized in Table [TJ Here in the static forcing model, our M2 tidal constituent has a period of half a 
day, and thus varies quite slowly with respect to the timestep {e.g. At = 0.5 s implies 86400 timesteps 
over T = 0.5 days). Nevertheless, one can see that because the energy of the solution satisfies a periodic 
forcing as well, the fixed tolerance regimes are only able to capture ~ 3500 timesteps of variation over 
the 86400 timesteps. This is because, regardless of the settings for e that one uses, the "energy" of the 
solution is only within the tolerance setting for ~ 1750 timesteps on either side of the peak amplitude. 
This behavior is demonstrated in Tableland Figure [3} The dioristic algorithms on the other hand, do a 
qualitatively better job of spreading out the p variation over the entire time domain [0, T). Moreover, as 
seen in Table[3[ the Type II algorithms owe a substantial component of their stability (and the robustness 
feature shown in Table [2]) to the fact that the de-enrichment functional is independent of t w , and thus 
is able to "catch instabilities" — so to speak — before they form; and hence they attempt to keep the 
average p value below a "stable setting," which by the CFL condition is a function of At. 

As a general trend, and as previously discussed, one might say that the Type I fixed tolerance p- 
enrichment scheme tends to p max in time assuming that the energy of the system is convex. Similarly the 
Type II enrichment schemes both have a tendency towards p m m in time under similar assumptions. The 
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Figure 4: Here we show the relative concentration of l after T = 2 days for the linear p = 1 solution. On 
the left is the solution with the static forcing, and on the right is the wavepool forcing. It is clear that 
after only two days, the wavepool has rejected a large amount of the contaminant that would otherwise 
be in the bay. 



Type I dioristic scheme, on the other hand, seems to demonstrate the most static behavior as a function 
of time, and it seems likely this will remain the case assuming that the signature behavior of the energy 
variation in 0, x [0, T) is not significantly dampened, whether the total energy functional behavior of the 
system is convex or concave or some mix of the two. 

It is also worth mentioning that the p-enrichment scheme does not come without a computational 
price, and in the Type I setting when e is small (i.e. sensitive) in a solution with very large gradients, 
the polynomial order rapidly approaches p max > so that when p is particularly high the p-enrichment 
scheme itself also becomes more computationally expensive. Generally we must be able to balance the 
computational cost of the enrichment scheme with the cost of increasing the local degrees of freedom of 
the solution, such that increased accuracy can be achieved without a resource cost that outweighs its 
worth. 

Finally let us briefly discuss the physical model. The wavepool forcing is designed with the intent 
of forcing the contaminant out of the bay by way of exploiting the mechanical energy of a hydraulic 
wavepool generator. Thus, we can simply test to see if, after 2 days, the amount of contaminant in the 
bay is reduced when employing the hydraulic wavepool generator, versus setting the static condition given 
by a f orce = m. We show the results after two days in Figure |4j First we note that the total (integrated) 
contaminant present in the domain in the case of the hydraulic wavepool model is 28% that of the case 
with no wavepool generator included. However, in the bay itself, the value is less than 18% that of the 
case with no wavepool generator included, clearly showing substantial declinature of the contaminant in 
only two days. This raises an interesting physical observation, which seems to suggest that fairly simply 
directed mechanical advection may offer a viable method of both declinature, and/or active transport of 
constituents which are inert with respect to the aqueous saline environments. 
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Figure 5: The left shows the mesh of the Brazos estuary for the eutrophication model, overlaid on a map 
of a Freeport, TX satellite image compliments to Google Maps, ©2009 Google, Map Data ©2009 Tele 
Atlas. The right gives the same mesh in terms of bathymetry, b = b(x). 



4.2 Neutralizing estuary eutrophication 

It has been generally observed that dead zones in estuarine outlets in the Gulf of Mexico correlate with 
algal blooms in the surrounding coastal regions. These blooms are believed to lead to dead zones in the 
local marine ecology by way of large population explosions in microorganism density which substantially 
deplete the oxygen in the surrounding sea water, in turn leading to devastating loss of local marine 
life. These blooms have been further correlated to excess nutrient levels (e.g. inorganic esters) in the 
gulf watershed leading to what are known as hypereutrophic coastal regions. In fact, the subsequent 
eutrophication in the local coastal region is known to be caused in part by excess relative concentrations 
of both phosphates and nitrates that are present largely due to the runoff pollution of man-made fertilizer 
in upstream agricultural centers spanning watershed riverbeds. 

Here we consider a test problem designed towards developing a method of neutralizing the rate limiting 
constituent (i.e. we choose phosphates as discussed in [18], since nitrogen is often in such vast excess) at 
an estuarine outlet. That is, we take the Brazos River estuary at Freeport, TX (see Figure [5]) and couple 
a remediation reaction to the contaminant flow from §5.2. 

Consider the chemically active form of the shallow water equations, which in conservation form sat- 

lsflGS ' 

d t ((Lj) + V x (Ht jU ) - Hsrfj = 0, 
d t (Hu)+V x & + S- v^x(Hu) - g(V x h = 0, 

© = (Hu ® u + X -g(H 2 - h 2 )^j , (4.4) 
/ n f n b \ 

= E(i - v i) '■■f-- n ^ - hr n ' 



given initial data 



reSK \ i=l i=l 



Ci=o = Co, u t=0 = u , t 3 -,|t=o = tj,o for je{l,2}, 
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Figure 6: Here we show the p level of the p =1-7 solution to the full eutrophication problem (with 
whirlpool and chemistry both active) after 28 seconds, using the Type I fixed tolerance scheme, with 
c = 1, c = 0.5 and At = 1 second. The z direction is mapped with respect to the bathymetry. 

where our equation is the same as that in §4.1, except for the fact that there are now two transported 
reactive chemical constituents u = Lj(t,x) for j E {1,2}. These constituents obey the law of mass 
action &fj = ^j(i), while the source term S = S(t,x) has been modified to include a whirlpool mixer 
= &{t,x) (as defined explicitly below) and a Coriolis parameter C = fuH for a constant coriolis 
coefficient / = 2$ sin 9, where = 7.292 x 10~5 rad • s _1 and 9 = 9{x) is the approximate latitude of the 
node. We add these terms {& and C) to the bottom friction already present in S from §4.1. 

The chemical mass action j^- may also be viewed as a source term in the transport equation where 
we have neglected the usual Fickian diffusion to a first approximation, as discussed in the introduction. 
Here, fc/ ir ,&b,r E M + are the forward and backward reaction rate constants, and Uj r ,Uj r E Z are the 
corresponding constant stoichiometric coefficients given an elementary reaction r in the reaction space 
9\. See [HJ [THl E01 |35l [38] for more details on these basic equations. 

Then, as a model problem we consider the aqueous (sea water) remediation reaction (buffered locally in 
the whirlpool to a pH ~ 8) of hardened gypsum and hydrogen phosphate into mineralized hydroxyapatite 
and salt, as observed in |23j by way of the proposed mechanism: 

CaS0 4 • 2H 2 + RP0 2 r — ^ Na 2 S0 4 + DAp, (4.5) 

4 NaOH 

where 

DAp = Caio-x(HP0 4 )x(P0 4 )6-x(OH) 2 _ x • nH 2 
is calcium deficient hydroxyapatite, and X E (0, 1]. 
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Figure 7: Each shows the solution at T = 29 days. The upper left shows the p = 1 solution to the 
inert transport problem, with chemistry and whirlpool off. The upper right is the p = 1 solution with 
chemistry on, whirlpool off. The bottom left is the p = 1 with chemistry on and whirlpool on. The bottom 
right shows the p 6 {1, . . . , 4} Type II fixed tolerance p-enriched solution with c = c = 1, chemistry and 
whirlpool turned on. 



Now, the gypsum CaS04 • 2H2O is a common and cheap industrial byproduct that demonstrates 
good hydraulicity, while hydroxyapatite is a primary constituent of bone, is osteogenic, non-toxic and 
potentially easy to reclaim from the environment (if such is even desired or necessary). The neutralization 
reaction of the sulfuric acid with the sodium hydroxide is taken to be diffusion limited, while the whirlpool 
pH is further buffered in the presence of sea water. The reaction rate of (4.5) is imprecise as extrapolated 
from |71 [21]. As a consequence, since the estuary water temperature varies annually from 12°C-30°C, 
we take as a weak estimate on the second order reaction rate assuming idealized conditions that kf ~ 
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1 L • mol 1 d 1 (though clearly this can be easily adjusted). 



Being concerned only with the reactants (4.5) may be partially decoupled from (4.4) such that they 
satisfy the following system of ODE's: 

d t t\ = -kfLii 2 and d t i2 = —kf titz, 

where i\ = [CaS04 • 2H2O] and 1% = [HPO4 - ] are relative concentrations. Applying the discontinuous 
Galerkin reactor scheme from |38j over our discrete timestep At = t n+l — t n we then arrive with: 

= tf e - k f t2At and = fie- k f LlA K 



where it becomes easy to reclaim the relative concentration t r - +l in (4.4) at every timestep. For more 
details on this family of DG chemical kinetic models see \6b\ I38| . 

The whirlpool mixer & = ff(t, x) sets the frequency of the gypsum release based on the local relative 
concentration of the phosphate, and also sets the speed at which the whirlpool rotates. As such, it 
satisfies the following condition at timestep n+1: 



sn+1 



L 2 ,ff for (sc 0l < x < x bl ) A (ti > L lf ff), , . 

Uff for (x ai < x < x bl ) A (u n < ttiim), 



where the whirlpool velocity is not allowed to exceed some limiting rotation velocity u\\ m which we set to 
u iim = 0.1 m-s -1 , restricted to the region (x ai , x a2 , x bl , x b2 ) at the heart of the whirlpool in the mouth of 
the Brazos River as shown in Figure [5] The limiting relative concentrations are given by = 1 x 1CP 8 
and 12 e = 1 x 10~ 2 , while the whirlpool velocity satisfies: 



Uff = r 1 m(x - x h ), for r = \/Jx - x h ) 2 , 

such that the velocity increases towards the heart of the whirlpool x^ and the direction of the rotation 
is determined by the vector m = (mi, m 2 ), setting m\ = —0.05 m-s -1 and 1112 = 0.05 m-s -1 . 

The corresponding boundary conditions are first split into Dirichlet conditions over the domain 
boundary dVth := dfl gu if U dft e stuary (see Figure [5]) corresponding to the Brazos River estuary and 
the Gulf of Mexico, except for now a free boundary subdomain is included dQi an d{t, x) and set such 
that d£li an d{t, x) = dQi an( i corresponds to the free boundary layer at the Freeport shoreline taken with 
respect to a coastal wetting and drying treatment (see below). Then by construction the total bound- 
ary fift o = Qhfi(t,x) at any timestep is given as dQh,o = d^iand U dClh where the remaining boundary 
conditions are given by: 

Cb = (Cb,l)land U (Cb,2,)gulf U (Cm) estuary > 
l>b = { L b,l)land U ( L b,2)gulf U (**, 3) estuary, 
Ub = («6,n,l • n + u b,r,l ■ T)land U («6, n ,2 ' n + u b,r,2 " T~) gu lf 
U i u b,n,3 ' n + u b,T,3 " T )estuary 

First, the Brazos River estuary is treated using the following no slip hydrogen phosphate inlet condi- 
tions: 

Cb,z\dn ei = C|an e , , WM j3 |,9Q e . = ^estuary COs(iO estu aryt) + 0.05 m • S _ , 

Uhr,3\dn e . = m • s" 1 and L b _ 3 \ d n. = 0.01, 



where we approximate the tidal frequency by setting uiestuary — 

0.729 x 10 _4 rad • s _1 with vr the tidal 
offset such that the variable amplitude per [22] achieves a mean stream velocity of a es t U ary = 0.7 m ■ s , 
and varies such that u bn> i\gci e . £ [0.02,0.12] m • s _1 . Finally a phosphate inlet relative concentration 
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of 0.01 is taken along the estuary outflow, where here we are only interested in the relative behavior, 
such that up to a rescaling all concentrations are admissible and may be viewed as unitless (e.g. volume 
fractions). 

Next, the open ocean conditions on the Gulf of Mexico boundary are set to the following: 

(b,2\dn Ei = T, u b>n>2 \dn ei = m • s~\ u b)T>2 \an ei = m • s _1 and Lb^dn.^ = 0, 

where the tidal constituents function T = T[t, x) is a linear combinations of its arguments such that 
T = T{K i, Oi, Qi, M 2 , S%, N2, K 2 ), where the arguments are the usual diurnal and semidiurnal tidal 
constituents. More precisely, for each constituent x £ {Ki> 0±, Q\, M 2 , S2, N2, K 2 } takes the form \ = 
nude cos(ojudet — vj) for u the tidal frequency, w = w(x) the relative offset, and aude = <Hide{&) the 
corresponding amplitude of each partcular constituent, each of which determined using the Eastern Pacific 
Tidal Database discussed in [33]. The remaining variables are quenched by the open ocean boundary. 
Finally, the land/beach free boundary at Freeport, TX are taken to satisfy: 

C&,i|an e . = u bnl \ d Q = u-grg, W6.r,i|an e . = uw® an d L b,i\dn e . = 0) 

where corresponds to the free surface and velocity components are determined by the sophisticated 

wetting and drying treatment presented in [8], while the chemical constituents are treated as the inert 
constituent is in §4.1. 

The wetting and drying treatment is beyond the scope of this paper, and we refer the reader to [8] 
for more details. It should however be noted that here, elements experiencing "wetting and drying" 
are restricted to the p = 1 case, as the treatment relies on a re-weighting of the slopes of lines (see 
j8]) in the case of a linear basis. This algorithm would require a nontrivial reworking to higher order p 
to fully achieve a homogeneous p > 1 solution. As such, the p-enrichment scheme is particularly well 
suited for just such a model, where one can readily p-enrich on interior elements by adapting within 

P € {1,.. . ,Pmax}- 

The numerical experiments show interesting behavior here. That is, the p level shown in Figure [6] 
follows the dynamics in the domain, where here we have used the high stability Type II fixed tolerance 
p-enrichment scheme. We see that along the edge of the tide, at the river inlet and at the boundary of 
the whirlpool that the p level drops to p m i n while elsewhere it adapts to the local perturbation in the 
solution, introducing a fair amount of additional structure into the polynomial basis. 

In the phosphate neutralization study the p = 1 cases were run at At = 2 seconds, while the adapting 
cases were run at At = 1 second. Some of the results are shown in Figure [7} Generally, we find that when 



the chemical reaction is active via (4.5) the amount of phosphates present in the Freeport coastal regime 
is reduced to 10% that of the solution when the chemistry is turned off. Interestingly, the addition of the 
whirlpool mixer (by which here we simply mean the vortical velocity field) does not appreciably improve 
the extent of the reaction, though it does have the effect of localizing the constituents to the mouth of 
the estuary and keeping the tidal fluxes from washing the constituents up/down the coast. The adapting 
case in Figure [7] shows very nice behavior here, giving sharper profiles along the edges of the contaminant 
flow in comparison to the p = 1 case, and indicates that the model is worthy of further study as a possible 
route to neutralizing eutrophication from fertilizer contaminants in gulf estuarine flows. 



§5 Conclusion 

We have presented a family of p-enrichment schemes designed for highly generalizable, computationally 
efficient, nonlinear free boundary type problems. This family of schemes was categorized into two classes. 

The first class are the fixed tolerance schemes. We have found that these schemes are particularly 
well-suited for problems whose "energy" is tightly bounded from both above and below, since the fixed 
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global tolerance otherwise is highly localized in either space or time. Moreover, the Type I schemes are 
attractive to many model applications which desire to emphasize particular "high energy regions" of the 
domain for enrichment, while the Type II schemes are appealing to application models where stability 
and computational efficiency trump the former concerns. 

The second class are the dioristic enrichment schemes. These schemes address the most obvious 
drawback of the fixed tolerance schemes, in that they adapt locally in time to the global spatial bounds 
of the regularity estimator. These schemes are well-suited for contexts which are complementary to those 
of the fixed tolerance schemes: namely, contexts where the "energy" varies substantially in both x and t. 

We then studied a pair of model problems. The first, a contaminant declinature model, was used to 
study the p-convergence of the solution, the accuracy and robustness behavior of the enrichment schemes, 
and the effect of weak entropy boundary layers. The physics of the model suggests that imparting 
mechanical wave energy might provide a simple way of diverting inert contaminants near important 
coastal regions. The second model was an estuary eutrophication study of a reactive multicomponent 
flow regime, where fertilizer runoff into a dead zone was counteracted by way of a neutralization reaction 
at the mouth of the estuary, and suggests that such a preliminary model might be able to provide a 
remediation mechanism in these areas. 

Future work is largely aimed towards model verification by way of employing the p-enrichment schemes 
to efficiently improve the accuracy of large applications models, such as hurricane storm surge. We also 
note that a very beautiful extension of the coupled chemical model from §4.2 is presented in [2], which 
includes the additional parameters of volatilization and sorption rates of chemicals in the atmosphere 
and suspended sediment, respectively, with the water. Studying the effects of these parameters on the 
eutrophication model, for example, might provide a more realistic understanding, as would verifying the 



hydraulicity and ecological impact of (4.5) in the context of sea water. 
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